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Abstract. We have designed a simple multi-scale method that identifies turbulent motions in hydrodynamical grid simulations. 
The method does not assume any a priori coherence scale to distinguish laminar and turbulent flows. Instead, the local mean 
velocity field around each cell is reconstructed with a multi-scale filtering technique, yielding the maximum scale of turbulent 
eddies by means of iterations. The method is robust, fast, and easily applicable to any grid simulation. We present here the 
application of this technique to the study of spatial and spectral properties of turbulence in the intra-cluster medium, measuring 
turbulent diffusion and anisotropy of the turbulent velocity field for a variety of driving mechanisms: a) accretion of matter in 
galaxy clusters (simulated with ENZO); b) sloshing motions around cool-cores (simulated with FLASH); c) jet outflows from 
active galactic nuclei, AGN, (simulated with FLASH). The turbulent velocities driven by matter accretion in galaxy clusters 
are mostly tangential in the inner regions (inside the cluster virial radius) and isotropic in regions close to the virial radius. The 
same is found for turbulence excited by cool-core sloshing, while the jet outflowing from AGN drives mostly radial turbulence 
motions near its sonic point and beyond. Turbulence leads to a diffusivity in the range Dturb ~ 10^' - 10^^" cm^ s"' in the 
intra-cluster medium. On average, the energetically dominant mechanism of turbulence driving in the intra cluster medium is 
represented by accretion of matter and major mergers during cluster evolution. 
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1. Introduction 

On many scales, astrophysical fluids show signs of turbu- 
lence whose dynamical contribution may range from signifi- 
cant, as in the case of the intra cluster medium (ICM) (e.g. 
Norman & Bryan Il999t iDolag et al.1 120051 : ISubramanian et al. l 
2006h to dominant, as in the case of the interstel- 



lar rnedium (ISM) (e.g. [L arson 1981; Goldreich & Sridhi 



Il995t iPadoan & Nordlundl 12002: Mac Low & Klesseni200' 
Turbulence is a fundamental ph enomenon that provides 



viscosity in accr etion disks (e.g. [Brandenburg et al.l 11995 



Balbus & Hawievlil99&). that transports matter in stellar atmo- 
spheres (e.g. lCanuto & Mazzitellilll99lh and mixes high- an d 
low-metallicity ICM in cluster cores (e.g. lRebusco et alJ 20061) . 

The direct numerical simulations of turbulence need to fol- 
low the turbulent cascade over a wide range of length scales. 
Recently, this has become feasible, as hydrodynamical simu- 
lations can reach fairly wide dyn amic ranges (of ~ 2-3 or- 



ders of magnitude in scales, e.g. iJones et alJl201 ll and refer- 
ences therein). However, these high-resolution simulations do 
not resolve the length scale of physical turbulent dissipation. 
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and subgrid turbulent closures that incorporate the evolution 
and effect of turbulence on unresolved scales have been devel- 
oped (e.g. ISchmidt et al. 2006t Scannapieco & Briiggen 2008 : 



Maieretal.ll2009t) 



The analysis of turbulence simulations of realistic systems 
(e.g. galaxy clusters) requires the separation of bulk and tur- 
bulent flows, and a number of strategies have been proposed 
in the recent past. A simple method would be that of comput- 
ing the turbulent velocity field as the residual respect to the 
ICM velocity field, averaged oyer spherical shells from the 
cluster centre Jn orman & Brvanl l 119991: llapichino & Niemeven 
20081 iLauet al.ll2009l) . Alternatively, on can compute the aver- 



age velocity field of the ICM via 3-D interpolation, and con- 
sider as tu rbulent the veloci t y structure below the interpola- 
tion sc ale /Po lag et all [20051: Ivazza e"tal]|2006l 120091 |2C)1 la : 
Valda rnini 201 ij). Alternative approaches focus on the decom- 



position of solenoidal and r otational comp onents of the ve- 
locity field dRvu et al.ll2008t IZ hu et al.'l201 Q), or empl oy sub- 
grid modelling (Ma ier et alj20 09: lapichino et al.l2011 ). These 
methods a priori assume limiting length scales of turbulence, 
possibly leading to inconsistent results. For instance, for a sim- 
ilar cluster mass the estimated amount of turbulent pressure in 
the cluster core may range from ~ 0.2 percent of the total gas 
pressure using sub-grid modelling estimates (Maier et al.2009) 
to ~ 2 percent of the total gas pressure by filtering the veloc- 
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ity field with a radial average (lapichino & Niemeyer 2008), 
to ~ 2 - 5 perce nt by using a filtering scale of ^ 300 kpc 
dVazza et alJl2009l) . 



In this article, we propose a method that locally determines 
the velocity coherent scale and uses this to distinguish the lam- 
inar and the turbulent components of the velocity field. Thus, 
the method needs no a priori assumptions of the typical scales 
of the flow in the simulated volume (Sec|2]l. The performances 
of our method are tested with an idealised setup in Sec l2.1l 
We present the first results on the properties of turbulence 
stiiTed by major mergers and accretions (Sec l3.1b . gas slosh- 
ing in cool cores (Sec l3.2l) . and active galactic nuclei (AGN) 
outbursts (Sec l3.3b . Our conclusions are given in Section|4l in 
the appendix we give an example of our algorithm in IDL 7.0 
syntax. 

2. Multi-scale filtering 

Turbulent fluids are characterised by a hierarchy of scales, 
ranging from the injection or driving scale, Lo, down to the 
physical dissipation scale, /tdis s, which sets the minimum scale 
available to the motion (e.g. iLandau & Lifshitz Il966t IShord 
In incompressible turbulence, the flux of kinetic energy 
across spatial scales is constant and, if the flow is stationary 
and uniform, the spectral energy distributi on for scales smalle r 
than Lo is described by E(k)dk oc k^^l^dk ( Kolmogo rovl 1 94 1 ) . 

This power-law translates into a simple relation between 
the physical size of turbulent eddies, I, and their internal veloc- 
ity dispersion, ctv: 



(1) 



(e.g. Landau & Lifshitz 1966). In most available numerical 
schemes, however, the smallest scale available to the fluid mo- 
tions is much larger than the physical dissipation one, thus 
breaking the power-law behaviour at some numerical s cale 
(~ 4 - 8 cells in most grid codes, IPorter & Woodward! e.g. 



1994 . 

The scale that conta ins the maxi mum kinetic energy is 
the "integral scale" (e.g. IShordl2007h . and for homogeneous 
isotropic turbulence it is given by 



A, 



2crl Jo 



E{k) 



dk. 



(2) 



while the largest correlation scale in the fluid. A, is defined by 
the maximum of k ■ E{k). 

According to this picture, the flow structure is uncorrected 
for scales » A, and the average velocity within this scale tends 
to the average fluid velocity. Based on that, we designed a re- 
cursive method to compute the average value of velocity around 
a cell for increasingly larger scales, until numerical conver- 
gence is achieved. The local mean field computed in this way 
(averaging for < A) can be used to compute the turbulent ve- 
locity fluctuations inside this scale. 

In detail, our algorithm works as follows: 
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Fig. 2. Trend with number of iterations of the mean local ve- 
locity field and turbulent velocity field reconstructed by our 
method for eight random points extracted in the first test of 
FigHI The top panel shows how the mean velocity and the tur- 
bulent velocity of each points change as the number of itera- 
tions is increased (the iterations are stopped when e < 10~^; 
the lower panel shows the trend of the fractional increase of the 
turbulent velocity field with the number of iterations (Eq|5]l for 
the same points. Our fiducial threshold value to stop the itera- 
tions, e = 0.1, as well as e = 0.05 and e - 0.01 is shown for 
comparison. 

- at a given n-th iteration, the components of the local mean 
velocity field around each cell are calculated as 

Xiir < L„)vi ■ Wi 



v(L„) 



(3) 



where Wi is a weighting function (e.g. gas density or gas 
mass); 
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Fig. 1. Maps of absolute value of the total velocity field (left, in arbitrary code units), of the turbulent velocity field reconstructed 
with our method (centre) and of the input turbulent velocity field (right) for two tests with a different slope for the background 
velocity profile (a) and for the number of injected turbulent patches, (see Sec l2.1l for details). In the top row we assumed 
a = -0.5, = 10 and o-y/vtot = 0.3, while in the bottom row we assumed a - -0.5, = 40 and cTy/vtot - 3. 



we compute the local "turbulent" velocity field at each n-th 
iteration: 



5v(Ln) = V - v(Ln); 



(4) 



the local and the turbulent velocity field are computed at 
each n-th iteration for increasing values of until the rel- 
ative variation of turbulent local velocity between two iter- 
ations is below the given tolerance parameter, e: 



(5v(Ln) - ^v(Ln-l) 
5v(Ln-l) 



< e. 



(5) 



In our case, AL - L„ - Ln-i is bound to be the minimum 
available cell size. The parameter e is a small tolerance pa- 
rameter, whose value is tuned by testing (we usually adopt 
e < 0.1, see below). 

Once convergence of Eq.|5]is reached, we fix A = and 
v(Ln) = Va, and we compute the local turbulent velocity 
field of the cell as 



6v - V - Va. 



(6) 



This procedure is repeated separately for each veloc- 
ity component. In simulations using the piecewise parabolic 
method (PPM) we set a minimum radius of - 4Ax at the 
start of iterations, since smaller scales can be affected by nu- 
merical dissipation. 

Note that if we choose the gas density as the weighting 
function, Wi, in Eq. [3] the numerical noise potentially arising 



by having gas cells at lower resolution far away from the cell 
location is minimised. This makes the scheme also readily ap- 
plicable to SPH, after interpolation onto a regular mesh. 

The numerical noise produced near strong shocks in the 
simulated volume affects the coiTect measurement of A and 
Va. In this case, the convergence of Eq|5]is made slower be- 
cause two different pre-shock and post-shock velocities are av- 
eraged across the shock. Strong shocks are characterised by a 
highly skewed distribution of velocities across the shock, and 
therefore monitoring the skewness of each velocity component 
in the volume around each cell is an efficient way of identify- 
ing the contribution from shocks to the local estimate of va. 
In detail, prior to our analysis we measure the skewness of the 
velocity field (separately for each component) in volumes of 
A^s = 8^ cells around each cell in the domain 



5i = 



1 



A^s 



Ns 



(Vj - Va)^ 



(7) 



where <Tyj is the variance of the velocity component inside the 
A^s volume around each cell. At each iteration of Eq. |5] we 
measure an average skewness inside the radius of integration 
by volume-averaging the previously measured values of Si, 
Sn = TjiS i/Nceii{Ln), whcrc A^ceii(in) IS the number of cells 
within the integration radius, L^. If a shock enters the inte- 
gration volume, the average skewness around the cell becomes 
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rapidly very large, and the iterations are stopped to avoid strong 
contaminations from velocity jumps at shocks. We found that 
stopping the iterations when 5„ > es, with es = 1 provides re- 
liable results for our simulated ICM with PPM methods. Using 
as a reference the realistic case of high-resolution ENZO sim- 
ulated galaxy clusters (as in Sec j3.1l l, we verified that slighly 
different choices in the range es = 0.5 - 3, or in the number of 
cells used for the local estimate of skewness, A^s = 5"* - 15^, 
yield very similar results (with differences at the ~ percent level 
on the values of turbulent energy) in the final 3-D distributions 
of turbulent velocity field with this method. 

In principle, more com plex and accurate shock-detecting 
schemes can be used (e.g. Vazza et al.ll201 Ibi and references 
therein), usually employing other physical quantities, e.g. gas 
temperature, pressure, sound speed. However, in our algo- 
rithm we ideally aim at reconstructing the turbulent field us- 
ing only the geometrical information on the velocity field. This 
makes the application of our method to a variety of simula- 
tions straightforward, without requiring fine tuning of parame- 
ters. The only internal parameters needed to stop the iterations 
of our algorithm are set after our preliminary testing: one for 
the convergence of the mean local velocity field {e - 0.1) and a 
second one to remove the spurious contribution of shock waves 
(fk - 1)- In a nearly homogeneous gas density distribution, the 
weighting function in Eq.[3]can also be omitted, and therefore 
the only physical field needed is the velocity field. 

Since our algorithm tries to reconstruct the typical scale of 
the signal in each point in space, this method is conceptually 
similar to the wa velet decornpositi on analysis used in turbu- 
lence studies (e.g. IMuzv et alJl99ll) . However, in our approach 
we do not aim at decomposing the 3D flows in i ts spatial com - 
ponents, as in the multi-resolution analysis (e.g. Mallatl[l989l) . 
but uniquely to constrain the largest outer scale of turbulence 
around each cell. 

In the appendix, we reproduce the source code of the basic 
version of the multi-scale filtering technique, written in IDL 
syntax for 3D distributions. 

2.1. Tests in two dimensions 

We tested our procedure and our choice of convergence pa- 
rameters against idealized 2D setups, in which we constructed 
combinations of average background velocity field and patches 
of chaotic turbulent velocity fields for a 200^ grid. For the reg- 
ular large-scale velocity field, we set up a simple radial inflow, 
according to VR(r) = A + B ■ r", where we set A = B = 1 
(in arbitrary code units). We tested a - -0.1, a = -0.5 and 
a = - 1 . These profiles represent a very generalised version of 
the profiles of radial velocity found in simulations of the ICM 
dNorman & Brvanlll999l: iFaltenbacher et alJl2005b . We added 
patches of chaotic velocity field by generating an additional 
2D velocity field, Vinput, with a random extraction from an en- 
ergy spectrum obeying the E{k) oc k^^^^ law. In our fiducial 
model, we imposed a minimum wavenumber kin =10 for the 
turbulent velocity; we also tested the cases kin = 1 and ki„ - 5. 
Then we randomly selected circular regions with random 
centres and radii. For the areas inside the extracted circular 



regions, we added the turbulent field to the background field, 

Vtot = Vo + Vinput. 

As an example, we show in FiglUthe maps of the total ve- 
locity field created in this way (left), of the turbulent velocity 
field reconstructed with our algorithm (centre) and of the input 
turbulent of velocity field (right) for two of our tests. Figure 
|2] shows the convergence of the local mean velocity and the 
turbulent velocity with iteration time steps for a random cells 
extracted in the first test of Fig[T] Our algorithm on average re- 
quires 5-10 iterations to converge on 6\ for each cell, within 
the e = 0.1 tolerance. The second panel of Fig|2]shows the be- 
haviour of the fractional change of the turbulent velocity fluc- 
tuation (Eq|5]l as a function of the number of iterations. In the 
vast majority of cases, stopping the iterations when this frac- 
tional change is below e = 0.1 represents a very good approxi- 
mation to constrain the turbulent field around the cell location. 
If we let the iterations proceed until a fractional variation less 
than e = 0.01 is reached, the number of iterations increases but 
the final improvement on the turbulent velocity field is not sub- 
stantial. Also based on our tests in the more realistic case pre- 
sented in Sec l3.1l we suggest that our fiducial choice of e = O.I 
is the best compromise between a robust reconstruction of the 
turbulent field and the speed of the algorithm. 

In Fig [3] we show the distribution of velocities for several 
tests, comparing the results of our algorithm to the input lami- 
nar and turbulent velocity fields. 

Our method performs well in reconstructing the original 
distribution of turbulent velocities in most cases, with no strong 
dependence on the background radial profiles (tests a)-c)). 
Misidentification can happen when the turbulent patches are 
so numerous that they frequently "merge" into a larger pattern, 
as in test d). In this case, the algorithm requires more iterations 
and a larger scan region to converge, and the estimated local 
mean field can be biased in this case. The best morphological 
reconstruction of turbulent patches is obtained when the tur- 
bulent structures are well separated, and their typical internal 
velocity is significantly different from the local mean velocity. 

Sharp cusps in steep velocity profiles can be misidentified 
as a turbulent fluctuation, as shown in test e), where no addi- 
tional turbulent field, Vinput, is added to the regular radial pro- 
file. However, such sharp peaks are unlikely in realistic simu- 
lations of the ICM. Our tests show that the use of e = 0. 1 is 
generally the best compromise between the need of a fast con- 
vergence of Eq.|5]outside of turbulent structures, and the neces- 
sity that the algorithm must not misidentify regular large-scale 
gradients as turbulent fluctuations. 

A second limitation of our method is that it relies on the 
assumption that the typical scale length of the laminar flow 
is larger than the maximum size of turbulent "eddies" in the 
simulated volume. When the two scales are comparable there 
is an excess of correlation within A, due to purely bulk mo- 
tions, which could bias high the estimated turbulent velocity. 
Indeed, when we impose kin = 5 or kin = 1 (tests g)-h) ) as an 
outer scale to add turbulent motions in our tests, differences are 
found between the distribution of velocity reconstructed by our 
method and the coiTect one. 

In Fig m we compare the power spectra of the input velocity 
field and the result of our algorithm for the few representative 
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Fig. 3. Distribution of the velocity modulus for the numerical tests described in Sec l2.1l From top left to bottom right: a) back- 
ground velocity profile with o- = -1 and = 10 additional patches of turbulent field with o-y/vo - 0.3; b) same as a), but with 
a - -0.5; c) same as a), but with a - -0.1; d) same as b), but with = 50; e) only background velocity field with a = -0.5, no 
turbulent patches; f) as in b), but with turbulent patches everywhere; g) as in b), but assuming an outer scale of turbulent motions 
kin - 5; h) as in b), but assuming an outer scale of turbulent motions ^in = 1. In each panel, the black lines show the total input 
velocity field, the blue lines the input turbulent field, and the red lines the field reconstructed with our algorithm. 



cases of tests b), d), e), f), g) and h). In intermittent as well 
as uniform turbulence cases (b) and f)) the input and measured 
spectra of turbulence are very similar for all scales, and they 
perfectly match the spectral shape of the input turbulence at 
the smallest scale. For the case without input turbulent field 
(test e)), as discussed above there is some residual pattern of 
misidentified turbulence. However, these patterns contain very 
low kinetic energy (~ 10"^ - 10""^ of the total energy), and 
they can be regarded as the unavoidable level of "noise" in our 
method. As mentioned above, for turbulent modes with a scale 
similar to the large-scale correlation of the profile of laminar 
motions (tests g-h)) our method faces a limitation, and at tur- 
bulent power spectra at the smallest k are not fully captured. 

We conclude that the method is accurate enough to sepa- 
rate laminar and turbulent pattern of motions in configurations 
similar to the simulated ICM. In realistic situations, the major 
caveat to the use of our procedure is that it is difficult to fully 
detect the largest-scale turbulent modes if they have a size sim- 
ilar to the largest scale in the computational domain. In this 
case an accurate measure of the outer injection scales of tur- 
bulent modes is only approximate, and the power spectra mea- 
sured are in general an underestimate at low k. This problem 
arises only when the physical injection scale and the physical 
scale of the ordered field are of the same order of magnitude. 
We will show that this unfortunate condition does not occur in 
the interesting cases of cluster mergers, cool-core sloshing and 
AGN-jets, which will be explored in the remainder of the paper. 



3. Applications 

We applied the multi-scale filter method to three important 
dynamical processes in galaxy clusters: turbulent motions in 
cosmological simulation of galaxy clusters, excited by merg- 
ers and accretion (Sect B.ll i; turbulence in sloshing cool cores 
(SectiT2ll; turbulence injected by AGN outflows (Sectl33Tl. 
In each case, we derived the turbulent power spectra, the 
anisotropy of the turbulent velocity and the resulting turbulent 
diffusion. 



3.1. Turbulence from mergers in galaxy clusters 

The presence of turbulent motions on scales » kpc in the ICM 
is inferred from several observations. Measures of Faraday 
rotation suggest the presence of chaotic super-Alfvenic mo- 
tions in the ICM, possibly excited by merger events (e.g . 



EnBlin & Vogt 



Bonafede et al 



F 

2003UMurgia et al J 120041: iGuidetti et all 12008 
20101: IVacca et al.ll2010h VIoreover, pseudo- 



pressure maps of cluster cores derived from X-ray observations 
and the lack of resonant scattering effects in the X-ray spectra 
provide hints of turbulence in the ICM ( Sch uecker et al.ll2004 



Churazov et al ] l2004l: ISanders & Fabianll2012h . Important con- 
straints on the fraction of turbulent and thermal energy in the 
cores of clusters are also based on the broadening of the hues in 



the em itted X-ray spectra of cool-core clusters ( Sanders et al 



20101) . In the next few years the satellite Astro-H with its high 
spectral resolution will provide an important tool to observa- 
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Fig. 4. Power spectra of velocity field for tests b), d) e) and f), g) and h) of Fig[3] The black lines show the total input velocity 
field, the blue lines the input turbulent field, and the red lines the field reconstructed with our algorithm. Note that the low k cut- 
off in the turbulent spectra of test h) is caused by our procedure of extracting small patches of turbulent motions in the simulated 
volume, while the original turbulent velocity field has by construction an unbroken power-law spectrum from k = 1 on. 



tionally constrain the en ergy ratio of turbulence in the ICM of 
real galaxy clusters (e.g. Zhuravleva et al.lEoi 1 ). 



Present-day cosmological numerical simulations routinely 
find that a significant amount of pressure support (i.e. ~ 10-30 
percent of the total pressure inside 0.57?vir) in the ICM is caused 
by chaotic m otions, coritinuously excited by major and minor 



mergers (e.g. Ilapichindl201 It IVazzall201 It I Jones et al.ll201 1 



for recent reviews). 

Here we study turbulence in the ICM of relaxed, merg- 
ing and post-merger galaxy clusters at high res o lution, with 



the set of s i mulatio ns presented in IVazza et alj (l2010ah and 



Vazza et al. ( 201 Iw . These runs were produced with the cos- 
mological adaptive mesh refinement code ENZO 1.5 (e.g. 
Norman et alJl2007tlCo"nins et al.ll2O10l) . ENZO is currently de- 
veloped by the Laboratory for Computational Astrophysics at 
the University of California in San Diego (http://lca.ucsd.edu). 

In Fig|5]we show the massive galaxy cluster. El, during a 
major merger event (z 0.6). We show the total velocity field in 
the centre of mass frame, the turbulent field after applying our 
multi-scale filter and the turbulent field below the fixed filtering 
length of 300 kpc or 1000 kpc. 

Inside the cluster atmosphere, chaotic motions are well de- 
veloped and significantly volume filling, and similar patterns 



of turbulence are detected regardless of the adopted scale for 
the filtering (as long as the filtering scales is a few ~ 100 kpc). 
This follows from the fact that, usually, the velocity field of the 
ICM in clusters is tangled for scales < 1 Mpc. However, the 
agreement between methods using a fixed filtering scale and 
our method becomes less satisfactory approaching the outer 
cluster regions, because towards 7?vir large-scale infall motions 
become more frequent, and large patches of laminar infalling 
gas are found in correlations with large-scale filaments. In these 
cases, equally strong smooth and chaotic flows can be found at 
roughly the same distance from the cluster centre, and distin- 
guishing among them is quite difficult if using a fixed scale. 

We compare in Figure |6] the mass-weighted radial profiles 
of total velocity and several estimates of turbulence in the same 
volume: by assuming fixed filtering scales (300 kpc and 1000 
kpc) and with our multi-scale filtering algorithm. In the same 
figure, we additionally show the results of slightly different 
choices of e to stop the iterations in Eq|5](e = 0.3, = 0.05 and 
- 0.01). Both fixed filtering scales yield in general a slightly 
higher turbulent velocity at all radii compared to our method; 
the differences with respect to the filter= 300 kpc are on aver- 
age very small. On the other hand, the turbulent velocity esti- 
mated with our method is lower by a factor ~ 2 -5 , with respect 
to the total velocity field of this cluster, for R > 500 kpc h'^. 
The choice of different values of e to stop the iterations of our 
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Fig. 5. Two-dimensional maps of total gas velocity fields through a cosmological AMR simulation at z = 0.6. Top left: total gas 
velocity (in [km s"']); top right: turbulent velocity field captured by our new multi-scale filter; bottom left: turbulent velocity 
field after the removal of L > 300 kpc scales; bottom right: turbulent velocity field after the removal of L > 1000 kpc scales. The 
side of each panel is 8 Mpc h"' . 



algorithm does not have dramatic consequences in the recon- 
structed velocity field. However, differently from the our tests 
in Sec l2.1l in this more realistic situation the choice of very low 
values of e may also cause problems, because to pin down the 
fractional change of the turbulent velocity field around the cell, 
very a large volume is scanned in the iterations, which can be as 



large as the cluster itself. Our fiducial choice of e = 0. 1 ensures 
a reasonable compromise between accuracy in the iterations, 
and the need of avoiding contamination from well-separated 
regions (and shocks) in simulated galaxy clusters. 

For completeness, we also show in FigQa map of the scale 
A of the velocity field inside the same scale as for the same re- 
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Fig. 6. Mass-weighted profiles of velocity from the centre of 
the cluster in Fig|5] The different lines show the total velocity 
(solid black), velocity field below the fixed scale of 1000 kpc 
(blue) and 300 kpc (light blue), turbulent velocity field recon- 
structed by our algorithm with fiducial parameters (e = 0.1, in 
red). We additionally show in grey the results of our algorithm 
for different choices of e to stop the iterations in Eq|5] e - 0.3 
(dashed), e = 0.05 (dot-dashed) and e = 0.01 (long-dashed). 




Fig. 7. Two-dimensional slice showing the distribution of co- 
herence scales (in kpc) for a velocity field of the same region 
as in Fig. 121 
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Fig. 8. Average radial profile for kinetic to thermal energy in- 
side radial shells for four simulated clusters, for the total ve- 
locity field (dotted lines) and for the turbulent velocity field 
reconstructed with our method (solid lines). 
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Fig. 9. Average radial profile of the outer scale for turbulence, 
A(r), for four simulated galaxy clusters. 



gion of Fig|5] The value of A shown is the average between the 
three velocity components. Across most of the cluster volume, 
A fluctuates in the range ~ 100-300 kpc. This suggests that, on 
average, the use of a simp le fixed filtering scale in this range (as 
perfo rmed in the past bv IDolag et al.ll2005t IVazza et al.ll2006 . 
2009h is still a good approximation to study turbulent motions 
in the ICM, and that the energy profiles obtained with this tech- 
nique in the past are fairly consistent with these new and more 
elaborated ones. 



F. Vazza, E. Roediger & Briiggen: A multi-scale method to measure turbulence. 



9 



1 

I 



1 H3- merging | H5 - posi-merger | 


11 32 74 158 327 661 


1327 


2672 


5330 





Fig. 10. Top panel: map of gas density for a slice of 100 kpc h ' 
through the centre of the major merger cluster H5 (right col- 
umn) and of the merging cluster H3 (left column). The top row 
shows the projected average gas density (in [p/pa.b], where 
Pci,b is the critical baryon density), the bottom row shows the 
projected map of the turbulent diffusion for the same regions 
(in units of [cm^ s"']). 



3.1 .1 . Turbulent energy budget in clusters 

We applied our filtering procedure to four additional galaxy 
clusters of total final virial mass Mft ~ 3 ■ 1 0"*Mo h'', which 



we already studied in previous works dVazza et al.l 1201 Obi) . 
These systems have different dynamical states: we have two 
"post-merger" objects (H5 and H6, with a merger with a mass 
ratio higher than 1 /3 for z < 0.5), one merging cluster (H3) 
and one relaxed system (HI, without evidence of past or on- 
going major merger for z < 0.5). Figure |8] shows the average 
radial profiles of turbulent and total kinetic energy within shells 
for these clusters, normalised to the thermal energy within the 
same volume, similar to Fig|6] The turbulent energy of cells 
is computed as iituib = Ax-'pcr^/2, where o-y is the modulus 
of the 3D velocity field below the assumed spatial scale, com- 
puted in Eq|6] and Ax is the resolution of the cell. Using the 
total velocity field would overestimate the turbulent budget by 
a factor ~ 3 - 10 at all radii for all objects. The ratio between 
turbulent and thermal energy flattens with radius in the range 
~ 0.5 - 2Ryii-, with iitmb/^^therm ~ 0.1, whilc the kinetic energy 
increases continuously with radius and approaches the thermal 
energy budget at Ry. The trend of /Stmb/^^therm with the dynam- 
ical state of host clusters is quite regular inside ~ 7?vh (post- 
merger systems present a higher content of turbulent energy 
compared to relaxed systems, while merging systems stay in 
between), whereas for < 0.5^vir the trend becomes sensitive to 
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Fig. 11. Average radial profile of the anisotropy parameter (for 
the total velocity field, in dashed, or for the turbulent velocity 
field, in solid) for four simulated galaxy clusters. 

the timing of the merger event, and to shock heating episodes 
(which affect the thermal energy of the ICM). 

The radial turbulent energy weighted profiles of the maxi- 
mum coherence scale, A(r), for these same clusters are shown 
in FiglH The trend of this scale with cluster dynamical states 
suggest that on average the most perturbed systems host the 
largest turbulent patterns, with A ~ 250 kpc, while the average 
values are almost half of that are found for the relaxed system 
HI . The differences tend to be smaller at /?vir, where on average 
all systems present correlation scales of A(r) ~ 100 - 150 kpc. 

These results cast doubts on the usual assumption of an 
injection scale of turbulence at shocks of the order of the 
curvature radius of accretion shocks, ~ 0.5 - 1 (e.g. 
Cavaliere et al.ll201 ih . For our clusters this would indeed im- 
ply values of ~ 1-2 Mpc h ' for the outer scale of tur- 
bulent motions. A, one order of magnitude larger than what 
we measure here. Also, the energy budget at ~ /?vir would be 
overestimated by a factor ~ 5 - 10 by assuming such a high 
value of A. The reason for this significant difference is likely 
that, in the outskirts of simulated galaxy clusters, turbulence 
is mostly injected at the scale ~ 100 - 300 kpc, typical of the 
density/pressure inhomogeneities of the ICM, downstream of 
accretion shocks. 



3.1 .2. Turbulent diffusion in clusters 

Our algorithm offers a straightforward estimate of turbulent dif- 
fusion in the simulated ICM, that is of the order of: 



turb 



0.11 ■ Ai -cr^ 



(8) 



where Ai 9A/20 (see for instance Dennis & ChandranI 



20051) ). Together with other mechanisms in the ICM, such as 
stellar feedback or galactic winds, turbulent diffusion may have 
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Fig. 12. Energy spectra of the 3D velocity field for three galaxy 
clusters (El, HI and H5). The dotted lines show the spectra of 
the total velocity field, the solid lines show the spectra of the 
turbulent velocity field. The gray lines shows the slope ofa-1 
and or = 2/3 to guide the eye. 



an important role in transpor ting metals from act ive galax- 
ies to the ICM. For instance, Rebusco et al. ( 2006h analysed 
the diffusion coefficient needed to model the metallicity ob- 
served around several central galaxies in nearby clusters, and 
found most likely values of the turbulent diffusion in the range 
Aurb ~ 10^** - lO^'^cm^ s-i. 



''turb 

We show in Fig[TO]the projected map of volume-weighted 
turbulent diffusion for the post-merger system H5 and for the 
merging cluster H3 along with their projected density for a 
slice of depth 200 kpc. In the major merger cluster H5 the ICM 
is characterised by volume-filling turbulent motions, which at- 
tain a maximum of Dtm-b ~ 0.5 - 1 ■ 10^"cm^ s"' in localised 
patches of a few ~ 100 kpc in size, and ~ 10^'' cm^ s"' else- 
where. In the pre-merger system H3 most of the cluster vol- 
ume is characterised by lower values of turbulent diffusion, 
Aurb < 5 • lO^'^cm- s"', but a few localised patches attaining 
higher diffusion are found related to minor mergers along the 
direction of the large-scale infall of matter (E-W direction in 
the image). Values in the same range are measured in the virial 
volume of the other two clusters. 

These results agree with the estimates of turbulent dif- 
fusi on of previous stud ies of the simulated ICM using trac- 
ers (iVazza et al. I l2010bh and imply that the efficiency of par- 
ticle transport in these simulated clu sters is much higher than 
that c aused by thermal diffusion (e.g. lShtvkovskiv & Gilfanov 
20101) . 

In the following sections we compare the distribution of 
turbulent diffusion in these runs with that of cool-core sloshing 
in the Virgo cluster i3.2i and with that of AGN outflows ( 13.3b . 



3.1 .3. Anisotropy of turbulence in clusters 

Our method also offers a way to monitor the anisotropy of the 
turbulent velocity field within the cluster volume, as a function 
of radius. 



m = 1 - 



v,tan 



20-2 ,' 

v.rad 



(9) 



where the turbulent field is decomposed into its tangential 
and radial component from the cluster centre. Knowing the 
anisotropy of the velocity field of galaxy clusters simulated in 
"simple physics" cosmological simulations is important to pin- 
point the additional effects of MHD instabilities, which can po- 
tentially lead to a radial alig nment of B and of the velocity field 



in the cluster outskirts (e .g. lOuataertI 12008c iRuszkowski et al 



20TlllParrish ^EoT Ij). Without a reliable method of detect- 



ing turbulent motions, laminar contributions of large-scale bulk 
flows can significantly bias the estimate of radial turbulent mo- 
tions in the ICM. 

Figure [TT] shows the average profiles of /3{r) for the four 
clusters; as a comparison we also overplot the corresponding 
profiles of the anisotropy parameter that we would obtain from 
the total unfiltered velocity field. Close to 7?vii, we note the 
striking feature that while the total velocity field is preferen- 
tially radial (J3 ~ 0.5 - 1) , the turbulent velocity field is close 
to isotropic (J3 ~ 0). Inside the virial radius, both velocity fields 
become preferentially tangential. Also in this case, the turbu- 
lent velocity field of all clusters is in general closer to isotropic 
than the total velocity field. These trends follow from the fact 
that while the bulk motions of satellites are characterised by 
strong radial motions towards the cluster centre, ram pressure 
stripping and the hydro-dynamical interaction with the ICM in- 
ject chaotic motions in a more isotropic way inside ~ /?vir- 

3.1 .4. Power spectra of turbulence in clusters 

Finally, we compute the power spectra of the turbulent velocity 
field and compare them with the power spectrum of the unfil 



tered velocit y field studied in our previous works (IVazza et al 



2010b[ l201 lal) . In Fig [12] we show the spectra for three clusters 
(El at z=0.6, HI and H5 at z=0) with a zero-padding technique 
and employing an apodisation function to avoid spurious ef- 
fects at the edges of the domain dVazza et al ]l2010bl:1 IValdarninil 
201 ll) . The wavenumbers are normalised to the virial radius of 
each cluster, ko = In/Ryii-. To determine the injection scale, for 
each cluster we plot the spectral energy per mode, k ■ E(k) (a 
Kolmogorov spectrum would have a oc k^^^^ scaling here, see 
the dotted gray lines in the Figure). Unfiltered spectra (dotted 
lines) present power law spectra for k > I - 2ko, with slopes in 
the range ofa~2/3-l. For the post-merger cluster H3, there 
is some flattening on scales of the virial radius. These results 
qualitatively agree with power spectra rep orted in the fitera- 



ture, based on different numeri cal methods (iDolag et al.ll2005 
Rvu et al.ll2()08tlXu et allEoOOl) . 



The turbulent velocity spectra of all clusters show a peak at 
k/ko ~ 10 - 30, corresponding to spatial scales ~ 0.1 - 0.3^vii- 
(~ 150 - 500 kpc h ' for the range of masses we consider 
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Fig. 13. Gas sloshing in the VIRGO cluster, triggered by a minor merger (i Roediger et al. 201 ll) . The images are slices through 
the centre of the simulation box, showing the absolute value of the total (top panels) and turbulent (lower panels) velocity for the 
three different resolutions: 2 kpc (LR), 1 kpc (MR) and 0.5 kpc (HR). The colour coding is in [km s '], each image has sides 
~ 250 X 300 kpc. 



here). We cannot detect a single sharp scale responsible of the 
injection of turbulent energy in the cluster volume, but a quite 
broad range of scales, consistent with the patchy distribution of 
scales shown in Figs|2]and|9l For the post-merger system H5, 
the turbulent energy peaks at a larger spatial scale compared to 
the relaxed cluster of equal mass, HI, implying the presence 
of strong turbulent motions caused by the most recent merger 
event. At smaller spatial scales, klh) > 50 (< 300kpc "'), 
aU clusters show a power-law spectrum with a slope slightly 
steeper than the Koknogorov one. 

Our findings suggest that, despite the clear power-law be- 
haviour of the spectral energy distribution of the ICM veloc- 
ity field (which runs for almost 2 orders of magnitude in spa- 



tial scales), turbulent motions indeed dominate the cascade of 
energy only for k/ko > 30 (corresponding to ~ 0.3Ryn, or 
~ 0.5 - 1 Mpc h"' for these masses). The spectral behaviour of 
the 3-D velocity field of the ICM for scales > 0.3 -0.5/?vir is on 
the other hand dominated by the pattern of velocities driven by 
large-scale infall, whose kinetic energy is mostly characterised 
by a laminar pattern. 

3.2. Cool core sloshing 

Our next example concerns ICM sloshing and the resulting for- 
mation of cold fronts (CPs). These are discontinuities in X- 
ray brightness and temperature, where the brighter side is also 
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the cooler one. They come in tw o varieties (see also review 
bv iMarkevitch & Vikhlininll2007h : merger CFs with stronger 
temperature contrasts across the fronts are the contact discon- 
tinuities between the i ntra-cluster media of tw o merging clus- 
ters, and sloshing CFs dMarkevitch et al.l200li) are named after 
their most likely origin. For the latter, the idea is that a gas- 
free subcluster moved through a galaxy cluster and its ICM. 
During the pericentre passage, the combined gravitational and 
hydrodynamical interaction slightly offsets the ICM in the clus- 
ter core without disrupting it. After the subcluster has passed 
the central region and moves away, the offset ICM falls back 
towards the main clust er DM peak a nd starts to slosh inside 
the main potential well ( Ascasibar & Markevitchl l2006 ) . Thus, 
sloshing CFs are contact discontinuities between gases of dif- 
ferent entropy, originating from different cluster radii. Usually, 
the subcluster passes the main cluster core at some distance, 
it transfers angular momentum to the ICM, and the sloshing 
takes on a spiral-like appearance, and so do the resulting CFs, 
which are wrapped around the cluster core Q. Measuring the 
amount of small-scale turbulent motions in these simulations 
is important because the excitation of turbulence around slosh- 
ing cool cores has recently been proposed as a mechanism to 
power radio mini-halos via turbulent r e-acceleration of y ~ 10^ 
electrons in the magnetised ICM (e.g. Mazzotta & Giacintuccil 
5008 ; ZuHone et al. 201ia). 

Roediger et al. ulated ICM sloshing specifically 



in the Virgo cluster, using the AMR code FLASH (version 3.2 
Dubev et al The simulations were performed in 3D in a 

simulation box of size of 3 x 3.5 x 3Mpc"'. 

Besides the superimposed sloshing and rotational large- 
scale motions, hydrodynamical instabilities at the CFs intro- 
duce a certain amount of turbulence, which we analyse here for 
the fiducial case with a subcluster of mass 2 • lO'^Mo, scale ra- 
dius 100 kpc, and pericentre distance of 100 kpc. To check the 
dependence of the turbulent velocity field on the numerical res- 
olution, we compared three re-simulations of the Virgo cluster 
with increasing maximum resolution: Ax ^ 0.5 kpc (HR run), 
^ 1 kpc (MR) and ^ 2 kpc (LR). 

Figure [13] shows the trend with resolution of the total ve- 
locity field (top panels) and turbulence (bottom panels) for thin 
cuts through the middle of runs LR, MR, and HR at the same 
time step. To run our algorithm over the same number of cells, 
we sampled all outputs at the resolution of the LR run (2 kpc). 

The increase of resolution enables us to capture the for- 
mation of Kelvin-Helmholtz (KH) instabilities along the spi- 
ral arms in detail, and to separate them more efficiently from 
the large scale rotation. The coherent rotation is characterised 
by values of ~ 300 - 500 km s"' at large scale in all runs. 
At low resolution (2 kpc) no KH rolls are observed along 
the spiral arms of the sloshing core, while at intermediate (1 
kpc) and high resolution (0.5 kpc) the KH instabilities can 
form, and develop patches of turbulent velocities of up to 
~ 100 -200 km s-'. 
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' This type of CF is reported to be ubiquitous 

dMarkevitch & Vikhliniij l2007h . and high-resolution obse rva- 
tions are available for several clusters ("see lRoediger et al]l201lL and 
references therein). 



Fig. 14. Power spectra for the total velocity field (dot-dashed 
lines) and for the turbulent velocity field (solid lines) for the 
Virgo simulation at the three resolutions. The spectra are com- 
puted for a cubic region with the side of 250 kpc. The x-axis is 
in unit of 27r/L, where L is the size of the box. The features at 
A; > 50 only present in run LR are an artifact of coarse resolu- 
tion. 

3.2.1 . Power spectra of turbulence in sloshing cool 
cores 

In FiglT4]we present the 3D power spectrum of the total and 
of the turbulent velocity field for a region of (500kpc)^ in the 
three runs. A large-scale correlation in the total velocity is 
found at all resolutions, with an overall slope steeper than the 
Kolmogorov spectrum. This correlation is mostly due to the 
correlation imposed by the large-scale pattern of rotation after 
the crossing of the cluster satellite. 

The increase of resolution creates a bump in the turbulent 
velocity spectrum for A: > 30 (corresponding to scales < 20 
kpc); this bump exactly corresponds to the peak of turbulent 
energy reconstructed by our algorithm, suggesting that the in- 
crease of resolution causes a real development of turbulent mo- 
tions in this range of scales. 

3.2.2. Anisotropy of turbulence in sloshing cool cores 

Furthermore, we computed the average radial profiles of the 
anisotropy parameter for the total and the turbulent velocity 
field in these three runs (FigfTSll. However, in contrast to our 
previous results on clusters, the tangential motions of the slosh- 
ing core strongly dominate the total velocity field inside 100 
kpc, and become more isotropic approaching the innermost 
cluster core. The convergence on [i{r) is very good at all radii 
> 20 kpc. The turbulent field is also preferentially tangential 
inside 100 kpc, but no clear convergence with resolution is ob- 
served within the cluster core, likely in response to the growth 
of KH instabilities at small scales as resolution is increased. In 
the top panels of Fig[T6]we show the components of the turbu- 
lent velocity field for a slice of 1 kpc in run HR. The turbulent 
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Fig. 15. Radial average profiles of the anisotropy parameter, 
yS(r), for the Virgo simulation at three resolutions. The dashed 
lines are for the total velocity field, the solid lines are for the 
turbulent velocity field. 



motions are weU-confined within the innermost 100 kpc of the 
cluster, also when seen in projection, as shown in the bottom 
left panel. 

3.2.3. Turbulent diffusion in slosliing cool cores 

As for the case of clusters, we computed the turbulent diffusion 
reconstructed with our method (bottom right panel of FigfTSll. 
Values of Dturb ~ I - 2 ■ lO^'^cm^ s"' are found in a thin stripe 
along the spiral arm of the sloshing, and in the innermost 100 
kpc around the cluster centre. It is intriguing that despite the 
differences in resolution and driving mechanism in the simula- 
tions, this range of values is very similar to what we obtained 
for the case of galaxy clusters simulated with ENZO (Sec j3.1b . 

3.3. AGN outflows in cluster cores 

Powerful outflows from AGN are a viable mechanism to in- 
duce significant turbulent motions in the innermost region of 
galaxy clusters, contributing to the mixing of metals in the 
ICM, and to the lifting of cold and low-entropy materials to the 
outer cluster vo lume, thus reducing or quenching cluster cool 
ing flows (e.g . 'Ciotti & Ostriker' 1997; C hurazov eTal] 12001 



Briiggen et ah.2005: Rebusco et al..,2006.) . Only recently fully 



cosmological grid simulations have reached a sufficient dy- 
namical range to model the evolution and feedback of AGN 



flie ICM (e.g. Xi 


1 et al.ll2009; Tevssier et alJl201 1 ; Dubois et al. 


2011;Kimetal] 


2011). 



Recently, the outflows from AGNs in the multiphase ICM 
has been simulated in detail with FLASH simulations by 
lGasnarietal.1 and ICasnari et al.l J^ah . These au- 



thors reported typical values of ~ 100 - 300 km s for the 
turbulent velocity injected by the AGN, ruling the onset of non- 




projeclion of SOO kpc j 



Fig. 17. Maps of turbulent velocity module (in units of 
[km s ']) for the a slice of 1 kpc centred on the Hydra run 
(left) and for the volume-weighted projection across 500 kpc 
(right). Each image has sides 300 x SOOkpc. 
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Fig. 18. Average radial profile of the anisotropy parameters of 
turbulent motions in the Hydra A simulation. We show as a 
dashed line the profile of the total velocity field, and as a solid 
line the profile of the turbulent velocity field. 



linear instabilities and leading to the condensation of cold gas 
filaments. 

We analysed the output of the simulated AGN-driven out- 
flow in the Hydra A cluster, simulated with FLASH 3.2. The 
volume around the jet injection has a side length of 1 Mpc 
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Fig. 16. Top panels: slices of 1 kpc through the centre of the Virgo run HR, showing the radial component of the turbulent velocity 
field (left, in [km/s]) and the tangential one (right,in [km/s]). Bottom panels: volume-averaged maps of turbulent velocity field 
(im [km/s]) for a line of sight of 200 kpc in run HR (left panel) and maps of turbulent diffusion for a slab of 20 kpc (right, in 
units of [cm^ s"']). Each image has sides 400 x 500 kpc. 



and employs AMR to reach the maximum spatial resolution 
of 0.5 kpc per cell. The AGN jet is reconstructed by two cir- 
cular back-to-back inflow boundaries 12 resolution elements in 
diameter (2 kpc). The jet material is injected in opposite direc- 
tions (at a velocity of vjet = 3 ■ 10"* km s"'), with a total power 



ofWi 



jet 



lO'*^ ergs-1. 



At the epoch analysed here, the bulk 



velocity along the jet is ~ 1500 - 1800 km s"', and a power- 
ful M ~ 1.3 shock has been driven into the suiTounding ICM. 
To mimic the observed offset between the shock centre and the 
AGN, a smooth velocity field of ~ 670 km s ' has also been 
imposed to the simulation, as a potential flow around a sphere 
of 100 kpc radius directed towards (-1,1,0). F or more details 



of the simulation setup, we refer the reader to iBruggenetal 
(l2007h and lSimionescu et all (l2009t) . 



In Fig [17] we show maps of the module of the turbulent ve- 
locity field for a slice of depth 1 kpc (left) and for the volume- 
weighted projection along 500 kpc (right) for a region around 
the cluster centre. 

As in the Virgo runs, our method is very efficient in re- 
moving the large-scale (> 50 kpc) laminar component of the 
velocity field, and highlights the complex pattern of turbulent 
structures associated to the interaction of the jet with the ICM 
atmosphere. Even if the bulk velocity along the jet axis can be 
as high as ~ 2000 km/s, turbulence is on average injected by 
rolls of size ~ 10-20 kpc, via hydrodynamical instabilities, 
at the low velocity of ~ 200 - 300 km s"'. Weaker motions 
(~ 10-100 km s ') are also injected in the innermost 100 
kpc, perpendicular to the jet axis, following lateral expanding 



F. Vazza, E. Roediger & Briiggen: A multi-scale method to measure turbulence. 



15 



weak shock waves. These findings are consistent with the re - 
cent ones of iGaspari et al.l (l2012bl) and lGaspari et al.l (l2012ah . 
Some of the turbulent features of the outer jet structure may be 
dete cted in nearby AGN inside the galaxy by A stro-H or Athena 
(e.g. lHeinz et"aLll2010l:lMendvgral et alfcoilh . 



3.3.1 . Anisotropy of turbulence in AGN outflows 

FigHDshows the average radial trend of y6(r) for the total (solid 
line) and for the turbulent velocity field reconstructed with our 
filter (dot-dashed). The large-scale velocity field is dominated 
by radial motions outside of 100 kpc (associated with the ex- 
pansion of the jet and of the running shock at the boundaries of 
the domain), with strong features of radial motions inside 100 
kpc, connected to high-velocity "knots" along the jet. The tur- 
bulent motions are close to isotropic along the whole jet struc- 
ture, suggesting that instabilities are very efficient in distribut- 
ing the driving kinetic power of the jet in 3D. 

This implies that the dissipation of kinetic energy in the sur- 
rounding ICM is very isotropic, which aids the AGN's capabil- 
ity of heating the cool core of the cluster. However, the presence 
of an even weak magnetic field carried with the jet is expected 
to have a sizable eff'ect, for instance affecting the exchange of 
heat between the jet region and the su rrounding ICM due to t he 
local magnetic pressure gradient (e.g. O'Neill & Jonesll2010l) . 
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Fig. 19. Power spectra for the total velocity field (thin upper 
lines) and for the turbulent velocity field (lower thick lines) for 
the same region as in Fig[T7] The different colour-coding shows 
the spectra for the component of velocities along the three axes 
of the simulation (direction 'Y' is the propagation axis of the 
jet). 



3.4. Power spectra of turbulence in AGN outflows 

As before, we used the power spectrum of the turbulent ve- 
locity field to investigate the important physical scales of the 
turbulence in the simulation volume (FigfT9]l. We also decom- 
posed the spectra showing the contribution from the three ve- 
locity components separately. While at very large scales there 
is some small excess of power along the initial injection axis of 
the jet (Y), as well as the lack of large-scale motions along the 
Z-axis (owing the large-scale velocity field along (-1,1,0) im- 
posed by construction), for most of the scales in the spectrum 
the power of the three components is very similar. Turbulence 
peaks at very small scales {k ~ 50 - 100, corresponding to 
~ 5 - 10 kpc), but the resolution of this run appears to be insuf- 
ficient to measure the spectral shape of the turbulent cascade in 
a clear way. 

3.5. Turbulent diffusion for AGN outflows 

Figure |20] shows the projected map of average turbulent dif- 
fusion across 500 kpc in the Hydra A run, measured as in 
the previous sections. The highest values of turbulent diffu- 
sion on large scales are associated with the most prominent 
turbulent rolls in the range ~ 50 - 100 kpc from the centre, 
with Aurb ~ 1.5 - 2 • lO^^cm^ s"', while much lower values 
(~ lO^^cm^ s"') are found in the inner 50-100 kpc around the 
cluster core. Extended patterns of efficient diffusion are also 
associated to the motions perpendicular to the jet axis in the 
downstream region of laterally expanding shock waves excited 
in the AGN outflows. The values we measure for the Hydra A 
run are compatible with the turbulent diffusion necessary to ex- 



plain the gradient of metallicity in the innermost cluster regions 
(Rebusco et al. 2006; Roediger et al. 2007), provided that what 
we measure here is turbulent diffusion ~ 160 Myr after the ini- 
tial jet was launched in the simulation box. 



4. Conclusion 

We have presented and tested a simple and robust algorithm for 
extracting the turbulent velocity field from a generic 3D veloc- 
ity field in grid simulations. The algorithm is based on an iter- 
ative geometrical analysis of the velocity field, and makes no 
a priori assumption on a physical scale when filtering out lam- 
inar motions. It iteratively calculates the local mean velocity 
and the turbulent velocity for increasing volumes around each 
cell of the simulation, until convergence is reached, taking into 
account spurious steep gradients related to shock waves. 

Our tests in Sec 12. II show that our method performs well 
in reconstructing the morphology and spectral features of in- 
termittent patterns of subsonic and transonic turbulent motions 
embedded in a large-scale background field. The main limita- 
tions of the algorithm emerge in the presence of cuspy back- 
ground velocity profiles (where the change of slope in the pro- 
file may be filtered as a small-scale turbulent fluctuations, if 
the variations occur on ~ few cells), and if the outer scale of 
turbulence is very close to the typical scale of the laminar flow. 

We applied our method to cosmological simulations of 
galaxy clusters with ENZO 1.5 (Sect l3.1b . to sloshing cool-core 
clusters (Sect l3.2l and to AGN outflows (Sect l3.3l ) simulated 
with FLASH 3.2. 
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Fig. 20. Map of average turbulent diffusion (in [cm^ s"']) for 
the projection along 500 kpc in the Hydra A run. The size of 
the image is the same as in FigfTTl 



In cosmological simulations, we find that turbulent veloc- 
ities are slightly tangential in the inner regions and isotropic 
in regions close to the virial radius. The same is found for 
turbulence excited by cool-core sloshing, while a jet produces 
slightly radial turbulence and isotropic turbulence near its sonic 
point and beyond. 

Our method naturally provides a way to estimate the tur- 
bulent diffusion in simulations, Dtui-b °^ A ■ ctv. We show in 
Fig|2T]the direct comparison of the volume and mass distribu- 
tion of turbulent diffusion in the three cases. Each mechanism 
presents a particular shape of the distribution. Turbulent diffu- 
sion from cluster mergers in general has a "simple" distribution 
with a maximum at Dtuib ~ 10^'cm^ s"' (volume-weighted dis- 



tribution), with the tendency of post-merger systems to present 
tails of enhanced diffusion, up to several ~ 10^"cm^ s"' . In the 
sloshing cool core, we observe two maxima in the distribution: 



one at ~ 10 cm s and associated with the innermost turbu- 
lent region close to the cluster centre, and one with less efficient 
diffusion, ~ 10^^ - lO^'' cm s , associated with the KH rolls 
along the spiral arms of the sloshing ICM. While the first fea- 
ture is very stable against the change in resolution, the second 
one evolves with the increase of resolution because of the ef- 
fect of a more efficient separation of differential rotation and 
turbulent KH rolls in our algorithm, and also because of the 
real onset of KH instabilities at smaller scales in the simula- 
tion. The turbulent diffusion in the Hydra A jet is expected to 
be more time-dependent compared to the other two. The distri- 
bution ~ 160 Myr after the jet launching presents a more com- 
plex distribution, owing to different patches of fast diffusion in 
the jet-lCM regions of interactions. 



Overall, the maximum values of turbulent diffusion attained 
by these mechanisms in the simulated ICM fall in the range 
^'turb ~ 10^^ - 10"'"cm^ s"'. Merger and accretion episodes in 
the ICM provide a more volume-filling mechanism of turbulent 
diffusion, and turbulent diffusion about one order of magnitude 
faster than jets and cool-core sloshing. We note that the average 
values we measure are of the order of the upper limits derived 
with XMM-Newton observations of pseudo-pressure fluctua- 
tions in Coma (Aurb < 3 ■ lO^^cm^ s"', Schueckeret al. 2004), 
provided that our 3-D distributions of turbulent diffusion are 
usually patchy, and make a direct comparison non-trivial. 

A physical ingredient that can alter some of our findings 
is the magnetic field. While its inclusion is not expected to 
change the dynamics of turbulent motions drive n by large-scale 
mergers and accretion on » kpc scales (e.g . IXu et al ]|200i; 
Ruszkowski et aklboi 1 ; Bonafede et aklEoil ), local amplifica- 
tion of B in shear flows can suppress the growth of instabilities 
and mixing motions a long the spiral arms of sl oshing structures 
jZuHone et al.ll2011b.) and along AGN-jets jO'Neill & Jones 
201 oh 



We end by noting that to resolve the turbulence excited by 
cluster mergers, sloshing and AGN-jets in the same simulation 
and keeping the hydro-dynamic details presented in these runs, 
one would need to cover scales ranging from 7?vir ~ 3 Mpc 
down to the presumed scale of physical dissipation at ~ 0.1 
kpc. As an illustration of that, we show in Fig.|22]a compos- 
ite power-spectrum of the simulated ICM, obtained by stack- 
ing the power spectra of total velocity and of turbulence for 
the simulations analysed in this work. Since the masses of the 
clusters as well as the kinetic energy input differ, we rescaled 
the turbulent energy to the turbulent energy contained within 
~ (200kpc)^ (which is well-captured in all runs) in the re- 
laxed cluster HI. This plot illustrates what might be the power 
spectrum for a cluster of total mass » 3 • IO^'^Mq and radius 
1900 kpc, subject to sloshing event and a low-power jet of 

100 Myr (the power of the jet 



jet 



lO**"^ erg s ' in the last 



is estimated from the amount of the rescaling needed to match 
the spectrum of the Hydra run to that of the simulated cluster 
HI at 200 kpc). 

The obvious caveat is that because the different turbulent 
motions come from independent volumes, the interaction of 
turbulent modes along the turbulent cascade cannot be cap- 
tured in this way. However, it is interesting that when rescaled 
in this way, the total velocity field from all simulations sits on 
a large-scale power-law, with k ■ E(k) oc k^'^/^ (equivalent to 
E{k) ~ k^''^^), from ~ 3 - 5 Mpc to ~ 0.5 - 1 kpc. Because in 
all simulations the background density profile is close to a beta- 
model, we propose that this broad-spectrum unveils the typical 
structure of correlated velocity field imposed in the stratified 
ICM, mainly because of geometrical reasons. When turbulent 
motions are intermittently injected into the ICM, they would 
add their spectral power to this pre-existing power law, causing 
bumps in the spectrum. They would also cause a flattening in 
the "background" velocity spectrum along their turbulent cas- 
cade, producing a spectral slope close to E(k) ~ k^^^^ - k^^. 
In the near future, the efficient design of filtering methods to 
analyse this wide range of dynamical scales will be an im- 
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Fig. 21. Volume (left) and mass (right) distribution of turbulent diffusion for all runs studied in this paper. We show in red the 
turbulent diffusion from the AGN-jet of Hydra, in blue the distributions of turbulent diffusion from the sloshing core in Virgo 
(we plot with different line-styles the distributions at different resolution, as in Sec l3.2b . and in black the turbulent diffusion from 
cosmological clusters (the different line-styles are for each different object studied in Sec j3.1l while the shadowed region shows 
the uncertainty in the overall cluster sample). 



portant challenge for the theoretical understanding of turbulent 
motions in cosmological simulations of large-scale structures. 
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Appendix 

Below we give an example of our multi-scale algorithm 
written in IDL 7.0, designed to recursively analyse one velocity 
component of a 3-D field. 

The program takes as input a 3-D velocity field, vel (we 
assume here it is written in binary format), as a regular grid of 
linear dimension of n cells. The fiducial tolerance parameters 
eps=0 . 1 and epssk=l are used in the computation, and a ker- 
nel scale of nk=8 is preliminarly used to compute the average 
skewness of the velocity field around each cell. 

The final output of the code are 3-D distributions of the 
turbulent field (turb), of the outer scale (scale), and of the 
skewness of the velocity field (sk). 

This code makes intensive use of the intrinsic IDL func- 
tions smooth, convol, and where to reduce the usage of 3-D 
loops. We verified that this greatly speeds up the execution. 

In Fig|T|we show the benchmark tests of the code applied 
to HR runs of the Virgo cluster (Sec l3.2l ) and interpolated to 
four different grid resolutions of 64^, 128^, 256^ and 512^'. For 
a sufficiently high resolution, the scaling of our algorithm is 
very close to linear with respect to the number of cells anal- 
ysed. The execution time to perform the multi-scale filter anal- 
ysis of one component of velocity in a 5 12-' grid is ~ 45 min- 
utes, and ~ 3 minutes for a 256^ grid. However, the details of 
the performance and scaling may change from problem to prob- 
lem, depending on the intermittency of the 3D velocity field 
under analysis. 

The source code for our algorithm can also be downloaded 
at this URL: 

jhttp://www.ira.inaf.it/~vazza/papers/turbofilter.pro" 
Also, a sample 3-D file of velocity extracted from a 
cluster simulation is given as an example, at this URL: 
|http://www.ira.inaf.it/~vazza/papers/sample-velocity.dat| 

pro turbofilter 

n=256 linear size of the grid 
vel=fltarr(n,n,n) ;...3-D array of velocity 
openr,3, 'vel.dat' ;...the grid is read 
readu, 3 , vel 
close , 3 

; . . .needed parameters and thresholds 
r2=uint(n""8 . 5-1 . ) ;.. .upper limit for L 
rl=4 ; . . . .lower limit 

turbo=fltarr(n,n,n) .turbulent field 

scale=fltarr(n,n,n) scale of the flow 

sk=fltarr(n,n,n) skewness 

scaleC* , * , *)= 0. initialization 

turbo (* , * , *)=8 . ; . . . 

eps=8.1 ;.. .tolerance in Eq.5 

nk=8 ; . . .number of cells to compute skewness 

epssk=l. ;.. .tolerance for the skewness 

drr=l. radial step for Eq.5 

; . . . preliminary computation of the skewness 

meanv = smooth (vel,nk,/EDGE_TRUNCATE) 

sc = absCCvx-meanv)/vel) 



10000 



1000 r 



100 r 



10 r 




1 I I I I I 

lO'' 10" 10"" 10" 10" 

NVd [cells] 

Fig. .1. Scaling between the CPU time employed by our algo- 
rithm and the total number of cells in the computational vol- 
ume for four different resolution of the Virgo HR run (Sec l3.2l i. 
Tests run on an Intel Quad-Core Xeon E5345 Linux Cluster 

kernel=MKE_ARRAY(nk,nk,nk, /float, value - 

1.) 

kernel (8, *,") = 0. 
kernel (*,0,") = 0. 
kernelC" , " ,8) = 0. 
kernel(nk-l,'' , '0 = 8. 
kernel ('Snk-l,-') = 8. 
kernel ("","nk-l) = 8. 

sk=convol C (vel -meanv) " 2 , kernel , /edge.truncate) 

sk=meanv"3/float(sk"l . 5) ; skewness, Eq.7 

scl = ; 

iterations to constrain turbulence 

for r=rl,r2,drr do begin 

width = 2."'r+l ; width of box 

meanv = smoothCvel .width, /EDGE.TRUNCATE) 
; . . .mean local velocity at each scale 

sc = abs((vel-meanv)/vel) ;.. .differential 
change in vel, Eq.5 

skm=smooth(sk , width , /EDGE.TRUNCATE) ; average 
skewness within L 

check of which cells are converged 

ibox=whereCCabs(sc-scl)/float(scl) It eps or 
abs(skm) gt epssk) and scalex eq 8.,nn) 

if nn gt 8 then begin 

turbo Cibox)=vel(ibox) -meanv (ibox) 
; . . .turbulent velocity in the cell 

scaleCibox) = float(r+8.81) ; . . .outer scale 

L 

endif 

scl=sc 

; . . zc=n"8 . 5 

;tvscl, [vel C zc) ,meanv(*,*,zc) , turbo C*,*,zc)] 
endfor 

; . . . .saves our final results 
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; . . .final turbulent field 
openw , 3 , " turb . dat " 
writeu,3,turb 
close, 3 

; . . .outer scale 

openw, 3 , "scale.dat" 

writeu, 3 , scale 

close , 3 

; . . . skewness 

openw, 3 , "skewness.dat" 

writeu, 3 , sk 

close, 3 

end 



